Introduction to epiAneufinder

Input formats

In the example above, a fragment file is supplied, but different input formats are possible for EpiAneufinder. Alternatives are a directory with bam files per cell or an already generated count matrix. The file type is identified automatically. The fragment file needs to be a single file with the ending tsv or bed. For the count matrix, a directory with three files is required, called matrix.mtx, barcodes.tsv and peaks.bed. The bam files need all to be in one directory with a bam index file each.

The filtering options are dependent on the input file type. In case, a fragment file is provided the parameter minFrags filters cells with too less fragments. In case, a directory with bam files is provided the parameter mapqFilter filters low quality reads. For the option to provide directly a count matrix, no filtering is currently implemented.

Note of caution We compared the results of CNV calling with the peak matrix compared to fragment file as input. We used the SNU601 dataset as in the publication and compared with the scWGS dataset as groundtruth. The correlation between the CNVs is worse when using the matrix as input (0.74 from 0.85). We loose reads that are not in peak regions so overall we get less coverage over the genome. The result is that we can miss some gains and losses that would otherwise have been identified with the fragment input. So we would suggest a bit of caution in the interpretation of the results, but the comparison with the groundtruth the method still gives valid results.

GC correction

The counts per windows are biased by the GC content of the associated genomic region, so that a GC correction step is essential. As this correction can be come time-intensive, we offer the user several options to balance runtime and accuracy in the parameter gc_correction.

The original method proposed in our publication, called loess is very accurate, but doesn’t scale well with large datasets. A faster, however more approximate alternative is the method called quadratic. Briefly, both versions fit a correction factor to adjust for GC-dependent count deviations, but using two different models. The option loess applies a LOESS model, while quadratic bins the data based on the GC content and fits a bin-wise quadratic polynomial model. Both versions are applied separately on each cell. Detailed descriptions and formulas can be found for loess in the epiAneufinder publication (Ramakrishnan et al., NatComms 2023) and for quadratic in the aneuFinder publication (Bakker et al, Genome Biology 2016).

A third alternative is to apply the LOESS function not separately on each cell, but once fit a model on the pseudobulk aggregates across all cells (using the median counts) and use this as a general correction factor for each cell. This option is called bulk_loess.

In evaluations on different datasets, we validated that the overall method performance differs only slightly depending on the chosen GC method. loess and bulk_loess were slightly more accurate overall. However, the option quadratic is the fastest, followed by bulk_loess. We would recommend for large datasets to apply quadratic (default option) or bulk_loess and for small datasets loess.

Other parameters

The example is run on the genome hg38, but any other genome available at BSgenome can be used correspondingly, after installing the matching bioconductor package. Corresponding to the genome version, the correct blacklist file need to be provided to the function and the corrected naming of the excluded chromosomes (variable exclude).

The window size parameter defines the region size for the CNV predictions, counts within each bin will be aggregated and one CNV annotation for each bin is obtained.

The minsizeCNV parameter regulates the minimal number of consecutive bins required for a CNV annotation, i.e. each CNV will contain at least minsizeCNV bins. Setting the variable to 0 disables the restriction that a CNV needs to have a certain size and also CNVs with one bin can be called.

The k parameters regulates the number of iterations of the segmentation algorithm. A maximum of 2^k CNVs can be found per chromsomes in k iterations. A large k parameter increases however the runtime substantially, so we set the default value to 4.

Split karyogram into subclones

Subclones in the dataset can be identified based on hierarchical clustering of the identified CNVs. The dendogram is cut at a chosen depth tree_depth.

#Load result table
res_table<-read.table("epiAneufinder_results/epiAneufinder_results/results_table.tsv")

#Need to reformat column names as dash are replaced by points in R
colnames(res_table)<-gsub("\\.","-",colnames(res_table))

subclones<-split_subclones(res_table,tree_depth=4,
                           plot_tree=TRUE,
                           plot_path="epiAneufinder_results/epiAneufinder_results/subclones.pdf",
                           plot_width=4,
                           plot_height=3)

#Data frame with the subclone split
head(subclones)
#>                                            cell subclone
#> cell-TCAGCTCTCGCGTTCT-1 cell-TCAGCTCTCGCGTTCT-1        1
#> cell-TACGCAAAGTCGTGAG-1 cell-TACGCAAAGTCGTGAG-1        2
#> cell-TAACTTCTCCCGGGTA-1 cell-TAACTTCTCCCGGGTA-1        3
#> cell-ACCATCCGTTCTGAAC-1 cell-ACCATCCGTTCTGAAC-1        3
#> cell-CTCAGCTAGAGACTCG-1 cell-CTCAGCTAGAGACTCG-1        3
#> cell-ACAATCGTCTCATCCG-1 cell-ACAATCGTCTCATCCG-1        4

Distribution into subclones

Annotated karyograms

Additional annotations can be added to the karyograms generated by epiAneufinder, for example cell type annotations or the subclonal distribution identified by split_subclones.

We recommend to save the karyogram as a png, defined via the variable plot_path, as pdf versions can get quite large.

#Load result table
res_table<-read.table("epiAneufinder_results/epiAneufinder_results/results_table.tsv")

#Need to reformat column names as dash are replaced by points in R
colnames(res_table)<-gsub("\\.","-",colnames(res_table))

annot_dt<-data.frame(cell=subclones$cell,
                     annot=paste0("Clone",subclones$subclone))


plot_karyo_annotated(res_table=res_table,
                     plot_path="epiAneufinder_results/epiAneufinder_results/karyo_annotated.png",
                     annot_dt=annot_dt)

Result karyogram with annotation

Read profile for individual cells

The predictions of epiAneufinder on individual cells can be further explored by evaluating the count profile of each cell. The set of cells for detailed investigation can be either subsampled randomly by assigning an integer value to the parameter selected_cells, then this number of cells are sample. Or a vector of barcodes can be given directly to selected_cells.

outdir<-"epiAneufinder_results/epiAneufinder_results/"

#Version 1 - Selected a specific set of barcodes, e.g. based on subclones
cell_barcodes<-subclones$cell[subclones$subclone==1]
cell_barcodes
#> [1] "cell-TCAGCTCTCGCGTTCT-1"

plot_single_cell_profile(outdir,
                         threshold_blacklist_bins=0.85,
                         selected_cells=cell_barcodes,
                         plot_path=paste0(outdir,"individual_tracks_clone1.pdf"),
                         plot_width=25,
                         plot_height=15)
#> quartz_off_screen 
#>                 2

#Version 2 - Randomly subsampling a certain number of cells
plot_single_cell_profile(outdir,
                         threshold_blacklist_bins=0.85,
                         selected_cells=3,
                         plot_path=paste0(outdir,"individual_tracks_random_subsample.pdf"),
                         plot_width=25,
                         plot_height=15)
#> quartz_off_screen 
#>                 2

Example track of a single cell